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The dynamics of capillary filling in the presence of chemically coated heterogeneous boundaries 
is investigated, both theoretically and numerically. In particular, by mapping the equations of front 
motion onto the dynamics of a dissipative driven oscillator, an analytical criterion for front pinning 
is derived, under the condition of diluteness of the coating spots. The criterion is tested against two 
dimensional Lattice Boltzmann simulations, and found to provide satisfactory agreement as long as 
the width of the front interface remains much thinner than the typical heterogeneity scale of the 
chemical coating. 
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I. INTRODUCTION 

The physics of capillary filling, originated with the pi- 
oneering works of Washburn [if and Lucas [|J , has pro- 
vided a constant source of interesting problems in fluid 
dynamics 0, 0] • Recently, with the burgeoning growth of 
theoretical, experimental and numerical works on micro- 
and nano-fluidics, the problem attracted a renewed in- 
terest di). Capillary filling is a typical "contact line" 
problem, in which the subtle non-hydrodynamic effects 
taking place at the contact point between liquid-gas and 
solid phase allow the interface to move, pulled by cap- 
illary forces and contrasted by viscous forces. Usually, 
only the late asymptotic stage is studied, leading to the 
well-known Lucas- Washburn law [l|, which predicts the 
following relation for the position z(t) of the moving in- 
terface inside the capillary: 



\t)- Z \0)= lHcOs{6) t 

3/Lt 



(1) 



where 7 is the surface tension between liquid and gas, 
is the static contact angle, fj, is the liquid dynamic viscos- 
ity, H is the channel height and the factor 3 depends on 
the geometry of the channel. Here, we focus on a two di- 
mensional geometry given by two infinite parallel plates, 
separated by a distance H (see fig. [TJV). 

One of the practical problems associated with the 
Lucas- Washburn law is the monotonically vanishing 
speed (Eq. ([1]) implies z ~ i -1 / 2 ) of the liquid front in 
the late stage of the process. The inertial effects that 
characterize the initial stage of the filling process, decay 
exponentially in time, thus leading to a higher sensitivity 



of the front motion to chemical or geometrical defects de- 
posited along the walls. As a result, the front often pins, 
since even the smallest hydrophobic spot represents a po- 
tential barrier that cannot be overcome by a front moving 
at a vanishingly small speed. 

In this work we deal with "chemical" coating, modelled 
through a position dependent contact angle, 9{z). In par- 
ticular, we focus on a dilute regime, where defects (hy- 
drophobic regions) are sparsely deposited on the walls, 
with a typical inter-distance larger than the length scale 
associated with inertial effects in the flow. This approach 
is interesting from the practical standpoint, because, in 
actual experiments, hydrophobic coating minimizes im- 
purities adsorption at the wall. Therefore, it is of impor- 
tance to provide a quantitative relation between the aver- 
age fraction and chemical properties of the hydrophobic 
surface, and the maximum length reachable by the inter- 
face along the capillary. Indeed, this immediately reflects 
into minimal requirements for the control of the coating 
process, or for the external supply needed to support 
front propagation up to a given distance. 



II. 



GEOMETRICAL SETUP 
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We consider a fluid penetrating a channel of length L 
with periodic hydrophobic spots, as depicted in Fig.[ljA 
As we shall see in the following, the periodicity is not 
an important restriction, as long as the distance between 
two consecutive hydrophobic spots is large enough, i.e. if 
the "diluteness" requirement is fulfilled. The channel is 
coated with an alternate sequence of strips A and B, of 
length wa and wb respectively, uniformly repeated along 
the channel length. Regions A and B have different wet- 
ting properties, characterized by the contact angles 9a 
and 8b- In the following, we consider the case 9a < tt/2 
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and Ob > t/2, a situation corresponding to a periodic 
sequence of attracting (hydrophylic) and repelling (hy- 
drophobic) sites (Fig. [33). 

In this Letter, we investigate the dependence of the 
front speed on the two lengths wa > wb and on the wet- 
ting angles Oa, Ob- In particular, we show that it is 
possible to reformulate the filling problem in terms of 
a damped-forced oscillator, and we develop a procedure 
that allows for an analytical prediction of the pinning po- 
sition. This analytical criterion is tested against direct 
numerical simulations based on the Lattice Boltzmann 
equation. 




FIG. 1: Cartoon of the heterogeneous coating. Top: both the 
walls of the channel are partitioned in an alternate sequence 
of strips A and B, with different chemical properties (8 A and 
Ob) and different lengths (wa and wb)- Bottom: force f(z) 
as a function of z. In the hydrophilic regions (A) the front 
is pushed rightwards, while in the hydrophobic ones (B) it is 
pulled leftwards. 



III. THE LUCAS- WASHBURN EQUATION 

The Lucas- Washburn equation describing a liquid of 
density pi and dynamic viscosity /i/, filling a capillary 
channel up to a distance z, reads as follows: 



d 2 ^ 
' dt 2 



-r\z 



dz 
~dt 



•/(*)• 



(2) 



The LHS represents the inertial terms, while the RHS 
reports the viscous and capillary forces. The parameter 
V = ~~fip m f ron t °f the viscous term defines a charac- 
teristic time, tfi = 77 -1 , beyond which the dynamics is 
controlled by the balance between the viscous damping 
and the capillary force f(z) — 2£££l^ll2_ Here H de- 
notes the height of the channel, supposed for simplicity 
2d. 

Eq. @ is exact only in the limit of liquid filling an 
empty capillary. In the case of a liquid-gas system, the 



density ratio of the two phases a = — must be taken 
into account [H, SH, H2| . 

As we shall see shortly, Eq. ([2]) stems from the balance 
between the total momentum change inside the capillary 
and the force acting on the liquid system. The explicit 
expression of the viscous drag is obtained assuming a 
developed Poiseuille profile inside the flow, with speed 
u(y) = 6jpy(H — y). In addition, the front speed is iden- 
tified with the average flow speed in the channel, namely 
z = u= l/H f" u(y)dy. 

For the case depicted in Fig. [TJ the forcing f(z) in the 
RHS of Eq. Q is a piecewise constant function, alter- 
nating between positive and negative values, as shown in 
Fig. [1)3. Upon introducing the new variable q = z 2 /2, 
Eq. ^ simplifies to: 



d q 
dt 2 



-v 



dq 
dt 



F(q) 



(3) 



where, due to the non-linear change of coordinates, F(q) 
is now a non-periodic piecewise constant function of q(z). 
Note that F(q) has the dimensions of a squared velocity, 
more precisely F(q) = 2V ca pV d cos(0(q)), where V cap = 
7 / ' hi is the capillary speed, and Vd — vi/H is the diffusive 
speed, being v\ = [iij pi the kinematic viscosity of the 
liquid. 

Eq. describes the motion of a forced damped os- 
cillator with total energy E(q) = K(q) + V(q), where 

K(q) = I fill an d V{q)> defined by the expression 

~llq ~ F(o)> are tne "kinetic" and "potential" energies, 
respectively. In the particular case F(q) = const (homo- 
geneous, non-coated channel), integrating Eq. © once, 
leads t: 

q(t) = Ft v + (q(t ) - Ft,) exp (~^^) ■ (4) 

which, in the limit t — to S> t v , tends to the asymptotic 
value 



q(oo) = Voo = Ft n 



"/Hcos(0) 
6pi 



(5) 



Taking into account the definition q = z 2 /2, it is readily 
checked that Eq. is equivalent to Eq. (TTJ. 



IV. CHEMICAL COATING 

Let us now focus the attention on the case of a di- 
lute dispersion of hydrophobic spots along the channel 
walls. In other words, we shall always assume that the 
length of all hydrophlic spots, wa is much larger than 
that of the hydrophobic ones, wb- Moreover, we also as- 
sume the natural condition that the length of hydrophilic 
spots is large enough to assure that inside each pulling 
region the front reaches its asymptotic Washburn-Lucas 
velocity, Uqo, given by expression ([5]). The critical mini- 
mal lenght, w A rlt , of the hydrophilic region which satisfies 
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the above requirement is easily estimated out of Eq. (JSJ . 
The requirement is that in a time At ~ 0(t v ) the front 
has travelled at least a distance Aq = (w^™*) 2 , leading 
to: 



(0) 



Typical values for water-air experiments are ^ = 10~ 3 
kg m s ,7 = 72 1CP 3 kg/s 2 , yielding a minimum 
sparseness requirement of the order of w c j[ lt ~ 10/im 
for H = 10/im. Of course, if the hydrophilic pulling re- 
gions are not long enough, the front dynamics is always 
in the inertial-transient region of the dynamical evolu- 
tion of Eq. f3|, and the overall dynamics becomes less 
universal and strongly dependent on the coating details. 



V. PINNING CRITERION 

Given the definition of q(z), we note that both the po- 
tential barrier to be overcome in order to "jump" over 
the hydrophobic obstacle, as well as the extension of the 
regions in which the front is accelerated (Fig. [2]), are in- 
creasing functions of z (and consequently of q). There- 
fore, if the front manages to approach its asymptotic ve- 
locity before encountering the first barrier, it will attain 
the same asymptotic velocity at the bottom of each of 
the subsequent potential wells. 

This observation opens the way to a simple predic- 
tion of the number of barriers the front is able to over- 
come before being pinned. Given that the kinetic energy 

— ^ v< ^' accumulated by the particle as it reaches 
the bottom of the wells, remains constant for all the n 
jumps, the final pinning position is fixed by the condi- 
tion < AV™ ax , where AV™ ax is the ever-increasing 
heigth of the n potential wells. 

Let us focus on the dynamics during the climbing 
of the n-th hydrophobic barrier, located between posi- 
tions (fcn-i and qi n (inset of Fig. [2]) and denoted by 
AV n (q) = V(q2n) — V(q). While the value of the poten- 
tial barrier decreases from AV™ ax — V{qin) — V{q2 n -i) 
to AV(q2n) — 0, kinetic energy is partly converted into 
potential energy and partly dissipated by viscous forces. 

If the coordinate q lies in the range qi n -\ < q < q2n, 
the temporal evolution of potential and kinetic energy 
obeys the following equations: 



dt 



AV n {q)=F B y/2K{q). 



(7) 



and 



-2T,K(q) + F B y/2K(q). (8) 

Equation © is obtained from Eq. by multiply- 
ing by q = y/2K(q). Dividing Eq. © by Eq. (0 and 
integrating, we obtain 



r-K 



dK(q) 



\/2K(l) JAV 



AV„ 



dAV n (q) 



(9) 
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FIG. 2: Potential V(q) as a function of q: the ratchet form 
is clearly visible, as well as the position-increasing potential 
barriers. 



where v!L 



6Mi 



This equation, expressing the 



kinetic energy if as a function of the potential barrier 
AV n (q), is easily solved, to yield 



AV n - AV 




(10) 
(11) 



where by definition = \J2K^. The front is pinned 
whenever two conditions are simultaneously met: kinetic 
energy is depleted, K — 0, and the front position lies 
inside the hydrophobic region, i.e. AV n > 0. 
This yields 



AV n = AVI 



log 1 



COS(£/b 



>0, 



(12) 

where AV™ ax = —Fs{q2n ~ <?2«-i) is the potential bar- 
rier of the n-th hydrophobic spot. Expressing AV£ 
an explicit function of n one obtains: 



as 



av: 



-F B w B 



n(wA + w B ) 



w B 



-F B w B L p 
(13) 

where L p = u(wa + wb) — wb is the pinning length. 
Inserting Eq. (TT3")) into Eq. (TT2")) . after some algebra, the 
following prediction for the dimensionlcss pinning length 
L p = L p / H is obtained: 



C a 



w B 



log 1 



cos(e A ) 
|cos(e B )| 



cos(8 A ) 

\cos(e B )\ 



(14) 



Here C c , 



Vcap / (72Vd) is a constant depending on 



the capillary parameters, w B = w B /H is the dimcn- 
sionless hydrophobic length, and —1 < cos(6 B ) < 0. 
Eq. p4[) uniquely identifies the pinning region as a func- 
tion of the coating properties and of the typical exten- 
sion of the hydrophobic spots, w B . As intuition sug- 
gests, the pinning length L p decreases for intense hy- 
drophobic coating, i.e. for large hydrophobic strength 
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(cos(9b) — > —1) and for high fraction of geometrical 
space covered by the hydrophobic component (large u>b)- 
Eq. (fT4|) highlights two distinguished limits, that is fast- 
fi Um 9 ( |cos(ef)| » °) and slow-filling ( rgg^ -> 0). In 
the former case Eq. (fl4)) delivers L p ~ C cap cos(8a)/wb, 
while in the latter the pinning length vanishes like L p ~ 
C cap cos(^) 2 /(2u;b |cos(#b)|). This asymptotic expres- 
sions show that in the fast filling regime the pinning 
length L p grows linearly with the hydrophilic cos(#a), 
independently of the hydrophobic strength, and stays fi- 
nite for any non-zero value of wb- On the other hand, in 
the slow-filling regime, the pinning length vanishes with 
the square of the hydrophilic cosf^) and is inversely pro- 
portional to the hydrophobic one cos(#s). Finally, in ei- 
ther cases, the pinning length scales inversely with the 
hydrophobic width wb- 

Eq. (TT4"|) also shows that the quantity 
L p wb/ (C C ap cos(6U)), is a universal function of 
cos(6a)/\ cos(6>b)| (see inset of FigJ3] ). Since this 
function is bounded in the range [0:1], Eq. fT5|) clearly 
shows that L p is always finite, unless wb = 0. 

Let us notice that in presence of heterogeneous di- 
lute coating (wa >> wb), one may use the previous 
expression as a first order guess of the pinning length 
Lp. For a typical case of a capillary filled with water 
at room temperature, of width H = lOfim, and with 
length of hydrophilic and hydrophobic spots of the order 
of w B = 0.01 - 0.1, we obtain L p ~ 50 - 500. 
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FIG. 3: Normalized pinning length L p as a function of wb 
for | cos(#a)/ cos(0b) ~ 10,20,30. The continuous line rep- 
resents the prediction given by Eq. (|16[1 . while the sym- 
bols are the result of the LB simulations. To match theory 
and simulation, we have taken Wg = wb — 8 in Eq. (|16|l . 
where 8 = 8 / H — 0.025. The other parameters used are: 
a = 0.028, L = 20, C cap = 5.7, 9 B = 92. Inset: Normalized 
pinning length L p as a function of | cos(0a)/cos(0b)|. Note 
the two distinct fast filling and slow filling regimes, as dis- 
cussed in the text. 
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TABLE I: Comparison between the prediction given by 
Eq. (|16|l . based on the generalised Washburn equation (|15|l . 
and the LB simulation, expressed in terms of the number of 
hydrophobic wells n cr i t = [L p /(wa + Wg) ~ w % ] that the 
front manages to overcome. The first column shows the ra- 
tio between the interface width 8 and channel height H. The 
simulations were performed taking into account the correction 
due to the finite width of the interface, = wb — 8. The 
second column displays the number of wells n^rit predicted 
by Eq. (|16jl , while the third one refers to the LB results nj!™ • 
Theory and simulations converge towards the same prediction 
in the limit 8/H < 0.03. 



VI. DISCUSSION 

The Lucas- Washburn equation ^ does not take into 
account two main sources of uncertainity. First, the inlet 
dynamics may be sensitive to the structure of the reser- 
voir, leading to significant deviations from the Poiseuille 
profile in the early stage of the filling process. Second, 
the description does not take into account the dynami- 
cal effects induced on the interface by the motion. On 
the other hand, it is well known that a moving inter- 
face may be significantly distorted by the viscous stress 
induced by the fluid motion, especially close to the con- 
tact line. Such a dynamical bending, may require the 
introduction of a dynamic contact angle (23l . [24], [26l \2l\ , 
leading to a driving force in the Washburn law, which 
is itself influenced by the interface motion, resulting in 
an unclosed "bootstrap" problem. Moreover, especially 
for hydrophilic coating, the shallow fluid wedge at the 
boundary may induce an extra dissipation [3], which 
modifies the dissipative terms. Asymptotically, the cap- 
illary speed becomes lower and lower, and one may argue 
that, at large times, the assumptions behind Eq.® are 
increasingly well fulfilled, thereby justifying the neglect 
of the above mentioned problems. Still, in many situ- 
ations the asymptotic "ideal" regime is never reached, 
and consequentely one needs to assess the influence of 
the previous effects on the filling process. 

To check the validity of the "pinning criterion" on a 
realistic system, we have resorted to a numerical sim- 
ulations of a two-phase fluid, using Lattice Boltzmann 
Eq uations [Tol [III [l3. Il3l. Ilij in two dimensions, to repro- 
duce the capillary filling in presence of coated patterns. 
The geometry is the same previously discussed (Fig. [!} , 
the only difference being the periodic boundary condi- 
tions imposed at the two lateral sides in order to ensure 
total conservation of mass inside the system. Similar 
problems, with or without heterogeneous coating, have 
also been studied recently in [28|, [23, H3, HH . 
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To consistently compare the simulation output with 
the theoretical prediction, we need to take into account 
the effect of two factors that differentiate simulations 
from the idealized description given by Eq. @: (i) the 
presence of a finite liquid-gas density ratio a = p g /pr, 
(ii) a finite width 8 of the liquid-gas interface in the sim- 
ulation. In order to take into account the unavoidable 
"resistance" of the gas occupying the capillary during the 
liquid invasion, one writes down the balance between the 
total momentum change inside the capillary and the force 
acting on the liquid+gas system. Simple calculations lead 
to a new equation governing the front dynamics [l9l l2l| : 



[aL + (l-a)z]^ 



+ (1-«)(S) 2 



(15) 



-r,[aL+(l-a)z]f t + f(z) 



where L is the total length of the channel and -q = 
^ H 5 (i+6Ai/g) ' being A; the slip length of the fluid . In 
the limit a —* and A; — > 0, one recovers the original 
Eq. (|2|). Equation (|15|) can be mapped onto the same 
Eq. ([3]), only with a slightly different definition of q = 
[aL+(l-a)z} 2 /2, and F(q) = 2(l-a)V cap V d cos(0 [*(?)]). 
The pinning criterion Eq. (I14p therefore simply changes 

to 







W B 



log 1 



cos(e A ) 

|cos(0 B )| 



cos(8 A ) 
\cos(8 B )\ 



a 



where C c , 



C cap (l + 6Vff) 2 /(l - a) 



(1-a) 

^ (16) 
and L is the 



dimensionless channel length L 



L/H. This generalised 
is tested against 



criterion, denoted from now on as L p h 
LB simulations. 

The above equation is exact if evaporation- 
condensation effects are negligible, i.e. when the 
gas is pushed out of the capillary, without any inter- 
action with the liquid. This is not the case for most 
mesoscopic models available in the literature [H, [l6j], 
which are based on a diffusive interface dynamics 
[32l [H, [13, [35| . As shown in a previous paper [2l[ , the 
dynamics given by (|15p is correctly recovered only in 
the limit of thin interface S/H — > and of negligible 
gas-liquid density ratio, a — > 

For the present case, the existence of a finite interface 
width 8 (generally around 5-6 grid points) in the simu- 
lation, has the effect of introducing an effective length 
of the hydrophobic obstacle wb- Therefore, in the fol- 
lowing, we compare the simulation results with the so- 
lution of Eq. (jTSJ) , with an effective length of the ob- 
stacle Wg — wb ~ S. Besides, in order to recover the 
correct hydrodynamic limit, the additional requirement 
S/H « 1 must be fullfilled (Table 1). 



VII. SIMULATION RESULTS 



Fig. [3J shows the prediction for the dimensionless pin- 
ning length Lu/ 1 as a function of the two main param- 
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FIG. 4: Time evolution of the front coordinate in the LB simu- 
lation (red points) and as numerical solution of Eq. (|15p (black 
line) . The simulations were performed using the following pa- 
rameters: a = 0.028, L = 3600,// = 200, wa = 144, w B = 
56,7 = 0.14,!/ = 1/6. For this case S/H = 0.025, in lattice 
units. The dashed line marks the prediction of Eq. (|16[) . In- 
set: comparison between imposed (black line) and observed 
(red dots) contact angles: note that while the imposed angle 
cos(#(z)) is a step-wise function, the simulation clearly shows 
an hysteresis [l7l. Il8j of the contact angle. The observed con- 
tact angle was obtained by interpolating the front shape with 
a circle of radius R, so that cos(#) = H/2R. Note that, due 
to the convex shape of the interface, the front coordinate (cal- 
culated at y = H/2) always precedes the location where the 
front first meets the hydrophobic region. 



eters wb and cos(£?a)/| cos(#b)| . The main figure dis- 
plays three sets of results, corresponding respectively to 
cos(9 a) / \ cos(9 b)\ = 10,20,30. The solid lines represent 
the numerical solutions of Eq. (|16|) . while the symbols 
correspond to the LB simulation results. As one can see, 
the numerical results are in satisfactory agreement with 
the theoretical prediction. To be noted that we have 
limited the numerical comparison up to pinning lenghts 
20 because of computational constraints. 



L max _ 15 



In the inset, the universal curve 



1 



is reported: it is possible to appreciate the two 
different regimes of slow filling, in which the 
product L p WB/(C cap cos(9ji)) grows linearly with 
cos(9a)/\ cos(#b)|, and fast filling, in which this quantity 
saturates to 1. 

It must be noted that the convergence towards L p h 
is accomplished only when the width S of the liquid- 
gas interface becomes thin enough with respect to the 
characteristic length H of the system (Table 1). In this 
regime, corresponding to S/H < 0.03, the overall asymp- 
totic trend and the final theoretical pinning position are 
in agreement with the simulation, as shown in Fig. 2] 
The discrepancy between the two curves during the ini- 
tial transient regime is due to the "vena contracta" phe- 
nomenon, reflecting the non-trivial matching between the 
reservoir and the capillary dynamics at the inlet of the 
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simulation setup [20| • Such small discrepancy can also be 
reabsorbed into an added mass term in LHS of Eq.0, 
as shown quantitatively in plj ]. 

At larger values of the ratio S/H, the error on L p be- 
comes higher, and in these cases the simulation system- 
atically gives L p > L p h (Table QJ . Such discrepancy is 
connected with the departure from hydrodynamics, for 
length scales of the order of the interface width, typical 
of any diffuse-interface description. 

The satisfactory agreement between the LB simula- 
tions and the hydrodynamical description based on the 
generalized Washburn equation (TT5|) indicates that ef- 
fects induced by the distortion of the interface profile 
are negigible, at least for the choice of coating properties 
discussed in this work. On the basis of simple hydrody- 
namical considerations, one would expect that the extra 
dissipation term induced by the wedge close to the con- 
tact line should introduce a correction p3l. [25l. l26j : 

F e *c<^logOff/A) (17) 

where A is the ultraviolet cutoff (inner length scale) up 
to which hydrodynamics applies. In LB simulations, it is 
reasonable to assume A ~ S. The correction term (fTT|). 
being localized at the interface, is always sub-leading at 
long times, differently from the viscous dissipation that 
is proportional to the volume occupied by the fluid. This 
extra-stress can change the local interface profile, in par- 
ticular close to the transition between two regions with 
different coating properties. This is shown in the inset 
of Fig. [U where the effective dynamic contact angle, is 
plotted versus the static (microscopic) contact angle im- 
posed by the boundary conditions. As seen, the dynamic 
profile is departing from the static one only when the 
front enters the hydrophobic region, i.e. where the in- 
terface tends to be pinned. Eq. (TTT)) suggests that extra 



dissipation has an effect also for long times, in the case 
of filling with very hydrophilic walls (6 a ~ 0°), and/or of 
filling in viscoelastic flows, depending on the rheological 
properties of the fluid. These interesting issues are left 
for future research. 



VIII. CONCLUSIONS 

In conclusion, we have derived an analytical criterion 
(Eq. (Unj)) to predict the pinning location of a moving 
front in a 2d plane geometry, based on the generalized 
Lucas- Washburn equation. This criterion has been de- 
rived under the assumption of dilutcness of hydrophobic 
spots, and has been satisfactorly tested against 2d Lattice 
Boltzmann simulations. We have also discussed limita- 
tions of both approaches: the first due to the neglect of 
interface distortion and extra-dissipation at the contact 
line for Washburn-like equations, and the second regard- 
ing the effects induced by deviation from hydrodynamics 
in LBE models. We have shown that, in the limit of fast 
filling process (6 a >> 0), and for sufficiently thin inter- 
faces, the analytical criterion and the LB simulations are 
in satisfactory agreement with each other. 

Work to extend the pinning criterion to more complex 
situations, such as in the presence of stochastic hetero- 
geneous coating (not diluted) and/or in the presence of 
geometrically coated interfaces, is currently underway. 
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